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We study the effect of viscoelastic dynamics on the frictional properties of a (mean held) spring- 
block system pulled on a rough surface by an external drive. When the drive moves at constant 
velocity V, two dynamical regimes are observed: at fast driving, above a critical threshold V^, the 
system slides at the drive velocity and displays a friction force with velocity weakening. Below 
the steady sliding becomes unstable and a stick-slip regime sets in. In the slide-hold-slide driving 
protocol, a peak of the friction force appears after the hold time and its amplitude increases with the 
hold duration. These observations are consistent with the frictional force encoded phenomenologi¬ 
cally in the rate-and-state equations. Our model gives a microscopical basis for such macroscopic 
description. 


I. INTRODUCTION 

Friction is behind many phenomena of our every day 
life experience such as the adhesion of car tires on the 
road [1, 2], the sound emitted by a violin [3, 4], or the 
wear of human articular joints [5]. Often, in these situ¬ 
ations, the two bodies in contact may display discontin¬ 
uous dynamics, the so-called stick-slip dynamics [6, 7], 
in which periods of rapid movement (slip) are followed 
by periods of relative rest (stick). At the macroscopic 
scale the stick-slip motion is justified by the empirical 
difference observed between the static and dynamic fric¬ 
tion coefficient, but a general explanation of this phe¬ 
nomenon from first principles is not yet available. The 
introduction of visco-elastic effects in the dynamics can 
be the key for a microscopic theory of friction. In fact, 
we have recently shown that viscoelastic solids sliding on 
a rough substrate can display a stick-slip instability [8]. 
Interestingly it was also pointed out that visco-elasticity 
is at the origin of the observed precursors, the micro-slips 
occurring at the onset of slip [9-11]. 

A scheme of the model studied in [8] is reproduced 
for convenience in Fig. 1. It consists of an ensemble of 
interacting blocks driven at velocity U on a rough sub¬ 
strate. The macroscopic friction arises from the real area 
of solid-substrate contact, which consists in the junctions 
between asperities. When the block is stuck on the rough 
substrate, the elastic energy of the spring /cq slowly accu¬ 
mulates over time, and is released when junctions break, 
letting the block move. The rupture of a single junction 
can trigger further rupture, with a characteristic time tq 
that characterizes the dynamics of the block-substrate 
contact. In the quasi-static limit (V —>■ 0+), when the 
driving time scale td is very slow compared to tq, the 
chain of events triggered by a single rupture can be very 



FIG. 1. (Color online) Sketch of the viscoelastic interface 
model, (a): the interface consists in blocks (labelled i,i + 
1,...) located at the positions hi, hi+i,... (empty squares) 
and bound together via a combination of springs (fci, /C 2 ) and 
dashpots {rju). Driving is performed via springs ko linked to 
the position w = Vt. (b,c): the asperities that provide the 
contact (highlighted with ellipses) hinder the sliding of the 
block on its substrate, with a random force fi'^{hi) acting as 
a microscopic static friction force, (c): the solid sled on the 
substrate and the asperities under stress are changed. 


large and is called an avalanche. 

In the case of purely elastic interactions between the 
blocks, the macroscopic friction force, namely the aver¬ 
age elongation of the spring ko, is the control parameter 
of the depinning transition, a second order transition be¬ 
tween a pinned phase and a phase where the system slides 
at the driving velocity V [12-16]. Visco-elastic interac¬ 
tions [17] change the nature of the moving phase, indue- 













2 


ing hysteretic behavior, as it was shown by Marchetti et 
aZ.[18, 19] in the context of the plastic depinning transi¬ 
tion of a vortex lattice. 

In the quasi-static driving limit (tq ^ td), the pres¬ 
ence of visco-elasticity (with a characteristic time t^) 
has been shown to produce a rich phenomenology of the 
avalanche dynamics, with main shocks, aftershocks and 
spatio-temporal correlations similar to that observed in 
seismology [8, 20-24], When the visco-elastic relaxation 
is very fast compared to the driving (r„ ^ td), these 
system-size avalanches can be understood as a stick-slip 
dynamics of the mean field model [8]. 

In the present paper we study the case where the driv¬ 
ing rate to is of the same order of magnitude as the 
viscoelastic time scale r^, but is still very slow compared 
to the duration of each single avalanche: tq ^ Tu — To- 
We show that in this limit, the system displays the ba¬ 
sic features of dry friction. In particular when a uni¬ 
form driving is applied, we observe a transition between 
a stick-slip regime (slow driving, <C td) and a steady- 
sliding regime (fast driving, 3> t^). When the driving 
velocity is not constant, as in the case of the slide-hold- 
slide protocol, our results display qualitatively the kind 
of behavior observed in the experiments [25], most re¬ 
markably, the existence of a peak in the friction force 
after the hold period, which becomes larger as the hold 
time increases. It must be noted that this kind of exper¬ 
iments on friction dynamics are usually modelled using 
the so-called rate-and-state equations [26, 27], which in¬ 
corporate the fact that the friction force does not simply 
depend on the relative velocity between the two nominal 
surfaces but also depends on the history of the contact, 
via a macroscopic “state” variable 9. In our model this 
dependence is encoded in the state of the microscopic 
visco-elastic elements. In this way our model gives a 
microscopical basis for the rate-and-state description of 
friction. 

The paper is organized as follows. In section II we 
briefly present the model. In section III we study the 
case of uniform driving, and complete this study with the 
slide-hold-slide protocol in section IV. We then interpret 
our results in terms of rate-and-state formalism in section 
V and conclude in section VI. 


II. MODEL 


The equations of motion for the model pictured in 
Fig. I are derived in app. A: 

r]odth^ = koiw - hi) + kiV'^hi + 

+ k2'V^hi — k2Ui (I) 

rjudtu^ = k2V‘^hi - k2Ui, ( 2 ) 

Here hi is the position of the block i, w = Vt is the 
driving term in the uniform driving protocol and tjq is 
the damping coefficient for the sliding of blocks rela- 



FIG. 2. (Color online) (a): Mean field model of Eqs. (1) and 
(2). (b): Equivalent model if the solution is stationary. 


tively to the substrate. The internal variable Ui ac¬ 
counts for the dashpots elongation, which induce the 
damping of the elastic force, — Ui), acting on 

the block i. The characteristic response time of blocks 
is To = 77o/max(fco, ^ 1 , ^ 2 )) while the characteristic re¬ 
adjustment time for the dashpots is Tu = r]u/k 2 - 

The random force f^'^{hi) mimics the contacts between 
the asperities of the block and those of the surface, and 
acts as a microscopic static friction force. When a contact 
is broken, two things happen: (i) the block moves and re¬ 
distributes stress to neighbouring blocks, (ii) the asperi¬ 
ties involved in the junction are renewed, as sketched in 
Fig. Ic. This renewal process makes the asperity land¬ 
scapes different for each block. For simplicity, we con¬ 
sider that the random force ff^^^ihi) acting on block i is 
completely independent from that acting on other blocks. 

Eqs. (I) and (2) were introduced in [8] and studied in 
d = 2 and in mean field for V = O'*'. Here we extend the 
mean field approximation (which corresponds to replace 
the laplacian with h — hi, where h is the average 

location of the blocks, see Fig. 2(a)) to the finite velocity 
case. Note that this model couples the Iii’s with h = 
hi/N, i.e. it actually represents N times the /c 2 , ??«, 
units. As it is usual in the mean field case, the values of 
the parameters are thus scaled as k 2 /N,f]u/N,ki/N, to 
ensure the extensive character of the total energy. We 
compute the macroscopic friction force 

a = ko{w — h) (3) 

accumulated in the system for different driving protocols 
and adopt the so-called “narrow wells” approximation 
[13]. In this scheme, the disorder ff'^{hi) is modelled as 
a collection of narrow pinning wells representing impuri¬ 
ties, with spacings z distributed as g{z) and with average 
z = Jq zg{z)dz. As the wells are very narrow, the disor¬ 
der force fi^ihi) that derives from this potential is 0 ev¬ 
erywhere except for countably many points. Within this 
approximation each block is pinned in a single narrow 
well (see the figure in app. A). For simplicity we consider 
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FIG. 3. (Color online) Evolution of the stress a over time for 
three velocities. At slow driving V = 0.001 (large amplitude), 

V = 5 (small amplitude) the stress oscillates periodically, 
with an amplitude that we denote Aa{V). At faster driving 

V = 9.8 (lower curve, almost always constant), the stress 
reaches a stationary value after a very short transient. The 
precise choice of initial condition only impacts the transient 
regimes. We used kg = 0.01, fci = 0.1, ^2 = 0.9 and Z = 0.1. 


narrow wells of constant depth (and shape), namely, the 
random thresholds are constant, = const. = 1. 

Our model provides a description of the interface be¬ 
tween two solid surfaces where all dynamical properties 
are concentrated in one of them, and the other (the sub¬ 
strate) is taken as inert. It has to be emphasized that 
the blocks have to be considered as microscopic single 
contacts between the surfaces, and as such they are in¬ 
dividually given a rather trivial and time-independent 
interaction law with the substrate. In addition, the ki, 
^ 2 , and rju elements (and to a certain extent also the fco 
springs) must be considered as part of the surface itself, 
the whole picture in Fig. 2 thus representing a small part 
of the two solid surfaces in contact. 

Note that here and in the rest of the paper, we give 
all quantities in dimensionless form. In order to restore 
physical units we need to re-introduce the units of the 
fundamental quantities: distance along the h direction 
(z), along the surface (x), force (f), and time (t), so for 
instance the previous unit value of means = If, 

spring constants ki,k 2 are given in units of fx^/z, fcg is 
in units of f/z, etc. 


III. UNIFORM DRIVING PROTOCOL 

At long times, under a driving performed at constant 
velocity V, the solution of Eqs. (1) and (2) becomes in¬ 
dependent of the initial condition and, in mean field, two 
phases are observed: for small V we hnd a limit cycle so¬ 
lution that corresponds to a stick-slip phase, whereas for 
large V a stationary solution exists. The stationary value 
of the friction force, cr(V), can be computed analytically 


for g{z) = S{z — z) using the crucial remark that if a 
steady sliding regime exists then in this regime h moves 
with velocity V. In particular, h = Vt — a{V)/kQ. In 
the stationary regime the /c 2 -plus-dashpot branch shown 
in Fig. 2(a) can also be thought of as connecting hi with 
w = Vt, instead of h, as indicated in Fig. 2(b). This is so 
because the additional stretching induced by the time in¬ 
dependent shift, w — h = a{V)/ko, is quickly absorbed by 
the dashpot, without altering the forces in any manner. 

The model in Fig. 2(b) coincides with the one studied 
by Dobrinevski, Le Doussal, and Wiese in [28] and can 
be solved exactly. In particular the friction force writes 
(see appendix C) 


a(V) = 



zk2 

gzk2/Vriu — X 


-(fco-l-fci), (4) 


and it is bounded by the two limiting values 

a(V ^ 0) = /*'“-I (fco + fci) (5) 

(t(U oo) = /*'"-I(fco + fci-h fe) (6) 

Note that : 

• The friction force decreases as the driving velocity 
increases, an effect called velocity weakening, and 
observed in tribology experiments for different ma¬ 
terials, especially at very low velocities [25, 29] 

• The velocity weakening displays a characteristic 
\/V decay at large V. 

• The dependence on fcp ^nd ki is limited to the last 
term in Eq. (4), that accounts for a constant shift 
of the whole cr{V) curve. 

Velocity weakening is a necessary link between the static 
{V = 0) and kinetic {V > 0) friction coefficients, and as 
such is known to be crucial in the triggering of instabil¬ 
ities in sliding systems, leading to stick-slip motion [6] 
and to the existence of earthquakes in sliding tectonic 
faults [30]. For a review on nanoscale models of friction 
and experimental results on nano-tribology one should 
consult [31] or the letter [32], which contains accessible 
references to the literature. In our model, velocity weak¬ 
ening is a direct consequence of the viscoelastic relax¬ 
ation. In fact, the model without viscoelasticity lacks 
any velocity dependence of the friction force, as in that 
case there is no internal time scale to compete with the 
driving time scale. It should be noted that the most com¬ 
monly observed velocity-weakening friction law is only 
logarithmic in V, but cannot be expected here because 
we introduced a single relaxation time scale, in a mean 
field model. This echoes the results found in [33], where 
a 1/U velocity weakening law is found to occur in a non 
interacting block model with random friction coefficients, 
where the relaxation time scale is present thanks to the 
non-instantaneous slips. One may also note that we do 
not expect any of the velocity-strengthening scenarios (as 
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FIG. 4. (Color online) Friction force against driving veloc¬ 
ity, for different values of feo (feo = 0.3,0.1,0.01,0.001,0.0001, 
from left to right); and constant fci = 0.1, k 2 = 0.9, ^ = 0.1. 
In the slow driving regime (on the left), motion occurs mostly 
through the abrupt slips during which stress falls from the 
stress Before slip, as (up), to the stress After slip, a a (down). 
In the faster driving regime, the dynamics is stationary and 
friction decreases with velocity. Crosses indicate the values 
of V for which the numerical integration has actually been 
performed. 


those proposed in [34]) to occur either, since we do not 
account for any of the faster mechanisms that become 
relevant under fast driving. 

We now turn on the numerical study of the mean field 
version of Eqs. (1) and (2) using an uncorrelated dis¬ 
tribution of pinning wells, with mean z, namely g{z) = 
z~^e~^/^. In practice we implement the Fokker-Planck 
method used in [8], adapted to the case of finite driving 
velocity. All technical details are left to Appendix B. 

In Fig. 3, the evolution of the stress over time is com¬ 
pared for three values of the velocity: at slow driving 
velocities the stress a{t) oscillates periodically with an 
amplitude denoted Act(F), while at fast driving a sta¬ 
tionary value is reached (i.e. A(t(F) = 0). Note that the 
amplitude Acr(y) of the oscillations is also the width of 
the stress drops or “gaps” that occur during the system- 
size events. 

This behavior points to a bifurcation of the dynamics 
as velocity is reduced. To study this effect, we report in 
Fig. 4 the maximum and minimum of the friction force 
over time, as a function of the driving velocity V, for var¬ 
ious values of kn- We observe that the stress gap vanishes 
smoothly at the transition point V = V^, pointing to a 
“second order” dynamical phase transition in the order 
parameter Act. In Fig. 5 we show the full phase diagram 
of the system in the k^-V plane. There, for various sets of 
(fci, ^ 2 ) we observe a divergence of as l/k^. This 1/fco 
dependence at low fcp comes directly from the uniform 
increase of the pulling force as k^Vt and implies that for 
any non zero values of ki,k 2 and V, there will always 



FIG. 5. (Color online) Phase diagram in the kg — V plane, 
for different sets of {k\,k 2 ) values: from bottom to top we 
used (fei = 0.0, fca = 1.0), (0.1,0.9), (0.5,0.5), (0.9, 0.1). In the 
lower-left region, the system behaves in a non-stationary way, 
i.e. we have stick-slip motion. For feo > fe 2 , the system is never 
non-stationary, even at quasistatic driving, as predicted in 
[35]. The bold dashed line has slope one, indicating a scaling 
of the critical velocity as ~ k^^ for small fee’s. The thin 
dashed vertical lines indicate the asymptotic behavior —>■ 0 
when feo —>■ fe^. 

be a value of feo below which some stick-slip dynamics 
occurs, even at very large velocities. 

In the stationary regime, our simulations confirm, over¬ 
all, the friction behavior of Fq. (4). In particular, in 
Fig. 6 we observe a clear velocity weakening with a char¬ 
acteristic 1 /V decay towards ct(F —>■ 00 ) = 1 — (feo + + 

k2)zl2. 


IV. SLIDE-HOLD-SLIDE PROTOCOL 

Slide-hold-slide experiments are an important tool in 
the investigation of the tribology of solids. When the 
sliding of a solid is interrupted for some time At, the 
contacts at the surface of the solid can strengthen over 
time, so that when sliding is resumed, a peak in the fric¬ 
tion force has to be overcome before one recovers the sta¬ 
tionary friction force. The amplitude of the friction peak 
increases with the hold time At since the relaxation is 
more effective when it has more time to act. Our model 
has all the necessary ingredients to reproduce the peak 
in the friction force after a hold period. In fact, when 
driven at a finite velocity, the viscoelastic elements do 
not have the time to completely relax to the most con¬ 
venient (i.e., lowest energy) configuration at each global 
position. If a hold time is given to the system, the me¬ 
chanical energy reduces as the viscoelastic elements relax. 
Upon resuming driving, this lower energy configuration 
requires a larger stress to initiate motion again. Thus the 
effect is expected to become stronger as the hold time in- 
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FIG. 6. (Color online) Velocity weakening in the steady slip 
regime (V > V^), far from the critical point. We observe that 
a decreases when the velocity increases. For U —>■ oo, the de¬ 
cay in the friction force towards the limiting large-velocity 
value (Jv=tx 3 goes as ~ V~^. Inset: plot in the log-log co¬ 
ordinates. The dashed line gives the pure power law with 
exponent —1. We used the same color code as in previous 
figures (curves that go to larger values are the lower fco’s). 


creases, saturating at hold times much larger than the 
viscoelastic relaxation time. 

To simulate this process, we set the driving velocity V 
to a value at which a steady sliding is observed. Once we 
reached a stationary dynamics, the driving is stopped for 
some time interval At and then resumed. The evolution 
of the friction force during this protocol in our model is 
shown in Fig. 7 (left). The height of the friction peak as 
a function of the hold time At is shown in Fig. 7 (right). 

We observe some similarities and some differences 
when comparing our results with experimental observa¬ 
tions. In experiments, during the hold time the stress 
shows a relaxation towards some lower value, something 
which is not observed in our model. The reason is 
that stress relaxation can occur only if some secondary 
avalanches (aftershocks) are triggered by the viscoelastic 
relaxation. When the pinning thresholds have a sin¬ 
gle value, then aftershocks are excluded and the decrease 
can not happen. However in models with a wide dis¬ 
tribution of pinning thresholds this relaxation has been 
observed [20] and is compatible with experiments. 

Qualitatively, the friction peak is similar in our simu¬ 
lations and in experiments: in both cases its height in¬ 
creases as a function of the hold time, but in experiments 
the increase is usually reported to be logarithmic, while 
here we observe an exponential saturation at large At. 
The reason for this difference is that our model contains 
only a single relaxation time constant, which naturally 
defines a typical time scale, above which the system is 
fully relaxed and thus no longer evolves. 


It has been long realized that the friction force cannot 
be described by a single valued function of the instan¬ 
taneous velocity. The history of the contact plays an 
important role in determining the actual friction force. 
Our model is an example of such a case, since the value 
of the friction force depends on the state of the viscoelas¬ 
tic elements, which in turn depend on the history of the 
system. We thus briefly recall the standard rate and state 
formalism, and then show how it provides an appropriate 
framework to understand our results. 


A. The RS Formalism 

Phenomenologically, the behavior of frictional contacts 
has been successfully described through the formalism 
named rate-and-state (RS) friction, originated in the 
works of Dieterich and Ruina [26, 27]. Instead of assum¬ 
ing that there exists only a kinetic and a static friction 
coefficient, the RS formalism assumes that the friction 
coefficient continuously depends on the relative velocity 
V between the two surfaces (the rate variable) and on a 
state variable usually called 0. We recall that by def¬ 
inition, the friction coefficient /r acts as a threshold for 
the friction force a actually arising from the contacts: we 
always have a < (j-Fn, where Fjv is the force normally 
applied on the solid. The usual RS form for ti{v,9) is: 

fj,{v,e) = tio + a\og{v) + b\og{e/Dc) (7) 

where a, b are positive constants. Note that we use a 
small V to indicate the instantaneous velocity, which may 
differ from the driving velocity V. The a term describes 
the so-called “direct effect”, the increase of friction with 
increase of the relative velocity commonly observed in 
many materials. The state variable 9 is supposed to fol¬ 
low a second equation, usually written: 


where Dc is the “critical slip distance”, i.e. the amount of 
slip (of the center of mass of the sliding block) necessary 
to break a newly formed junction. 

Under steady sliding, we have 9 = 0, 9 = Dc/v, and 
the RS equations (7), (8) simplify into: 

ti{v,9) = fio +{a-b)\og{v). (9) 

If & > a this equation describes the phenomenon of ve¬ 
locity weakening, namely a reduction of the friction co- 
efhcient when velocity increases. Velocity weakening is a 
crucial ingredient involved in the description of seismic 
phenomena [30]. In the case in which the contact is at rest 
{v = 0), we get 9 = 1, i.e., 9{t) = 9^ +t, and according 
to (7) we obtain an increase of the static friction coefh- 
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FIG. 7. (Color online) Slide-hold-slide experiment. Left: Evolution of the stress o over time under this protocol, around the 
time when driving is resumed {t = 0.015 here). The driving velocity is 17 = 1000. 

Right: For many different hold times At, we record the maximal deviation (A(j)maa: from stationary stress (blue points). This 
deviation corresponds to the gap between static and kinetic friction coefficient. In our model we observe a linear growth of the 
gap with At while in experiments a logarithmic increase is measured. The dashed line gives the purely exponential saturation: 
Arr max — A(joo(l — e with T — 2 and Actoo = A(jmax(At = oo). 


cient with the time of contact. RS equations have been 
used to describe the behaviors of frictional systems un¬ 
der a variety of non-steady sliding conditions, providing 
an excellent phenomenological description for numerous 
systems. 


B. Interpretation of our model in terms of RS 
Equations 

Under uniform driving, our model displays two features 
that are common to many frictional systems: a velocity 
weakening friction law, and a bifurcation to a stick-slip 
regime for low driving velocities. We have described in 
Section III how the velocity weakening effect appears in 
the model. Now we will see that assuming the existence 
of this effect, the bifurcation to a stick-slip regime at low 
velocities is captured in all its detail by the use of the RS 
formalism. 

The usually assumed form Eq. (9) of the friction coeffi¬ 
cient is not appropriate to accurately describe the results 
observed in our model in the steady sliding regime. On 
the one hand, our model does not have any ingredient 
that could produce a direct effect in the RS equation 
(7), suggesting that we should use a = 0. On the other 
hand, the logarithmic b term is not appropriate in our 
case, mainly because we have considered only a relax¬ 
ation mechanism that is described by a single time con¬ 
stant, the one associated to the dashpot and ^2 springs. 
Instead, from our simulations of Sec. Ill (as for instance 
those in the right part of Fig. 4) we can construct a sim¬ 
ple analytical expression for the friction coefficient in the 


steady sliding regime: 

^steady I (10) 

where the numerical values correspond to ki = 0.1, k 2 = 
0.9 and z = 0.1. We remark that this is just a simple, 
although accurate fitting to the numerical results that 
allows the analytical treatment presented below, but it 
has no additional particular significance. 

We now assume that under any non-steady sliding sit¬ 
uation, the value of the friction coefficient can still be 
written in terms of the function found in the steady state. 
Using the fact that in the steady state v = Dc/9, we 
write: 

= + + + ( 11 ) 

The time evolution of 9 will be assumed to be described 
by the standard RS evolution, Eq. (8). The sliding con¬ 
tact is coupled to a driving spring pulled at constant 
velocity V, so that the pulling force reads: 

tr = {V t — x)ko (12) 

where x is the spatial average coordinate of the contact, 
i.e. X = V. This force a is also the actual instanta¬ 
neous friction force arising from the contacts. Our aim is 
now to determine the temporal evolution of x and cr and 
reproduce the bifurcation behavior that we obtained in 
Section III. 

If steady sliding is assumed, the force must bal¬ 
ance the pulling force, and we get 9 = Dc/V, v = V. 
However, steady sliding may be unstable due to the 
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FIG. 8. Friction force as a function of velocity for the same 
set of values as in Fig. 4, using the RS formalism. The form 
of the steady sliding part of the curves is introduced by hand. 
We used a value Dc = 0.06 to achieve the best fitting. The 
RS formalism yields the form of a max and amin for V < V^. 


following process: if at some moment x is stopped, 
the pulling force starts to increase according to u = 
tT° + Vtko- On the other hand, the friction force starts 
to increase due to the increasing of 0, according to 
= fJ-{0o + t). If the initial increase of is more 
rapid than the increase of tr, the contact will stay stuck 
as long as > cr, until eventually the pulling force 

becomes larger than the pinning force, and a rapid slip 
moves the value of a: to a new position. If /rfjv < cr, slip 
occurs and is essentially instantaneous since we have not 
added inertia to the contact. During the slip stage, the 
value of 0 is reduced following Eq. (8), that in the rapid 
slip limit can be written as d0 — —0dx/Dc and analyti¬ 
cally integrated. Slip finishes when a drops below /iF/v 
again, and a new stick slip cycle begins. We have inte¬ 
grated Eqs. (11), (8), (12) and in fact obtained steady 
sliding for large V, and a stick slip behavior at small V, 
in which the value of cr oscillates between two values Umax 
and (Tmin- In Fig. 8 we plot the values of Umax and CTmin 
(or the single value a in case of steady sliding) as a func¬ 
tion of V for the same set of kg values used in Fig. 4. It 
can be seen that the value of Dc enters in the problem 
only in combination with the spring stiffness, as kDc- In 
Fig. 8 we used = 0.06 to achieve the best fitting 
with the results of Fig. 4. We note that as suggested by 
its original meaning, Fc is of the order of z, which is the 
distance at which the correlations between pinning forces 
disappears. In addition to the coincidence of the two fig¬ 
ures in the steady sliding regime (which is enforced by 
hand through the choice of the function in Eq. 

(10)), we see that the RS formalism gives a very good 
coincidence in the low velocity, stick-slip regime. 

A careful analysis of our RS equations near the bifur¬ 
cation point shows that the two values amax and amin 
depart symmetrically from the branch of steady sliding. 


and are such that amax — amin ~ {V‘^ — which 

corresponds to a Hopf bifurcation. Note however that 
the validity range of these scaling becomes progressively 
smaller as kg is reduced, and in the fco —t 0 limit we get 
the scaling amin ^ const, Umax - amin ~ (R^ - R). 

VI. CONCLUSION 

We have presented a detailed analytical and numerical 
analysis of a viscoelastic model of friction in mean field 
approximation, in which the driving velocity competes 
with the time scale of the viscoelastic effects within the 
system. 

Our main findings are the following. At low driving ve¬ 
locities, we obtain a stick-slip dynamics with amplitudes 
of the stress oscillations that decrease with increasing 
driving velocity {V) or increasing driving stiffness {kg). 
Beyond a certain critical driving velocity {V > U'^), the 
amplitude of these oscillations becomes zero, i.e. the slid¬ 
ing is smooth. The transition between stick-slip and 
smooth sliding occurs in a continuous manner and is 
very well described by a phenomenological rate-and-state 
analysis. In the smooth sliding regime the friction force 
reduces as a function of the driving velocity, reproduc¬ 
ing the well known velocity-weakening phenomenology. 
In our model, this effect is originated in the existence of 
viscoelastic elements that set an additional time scale for 
the dynamics. Finally, the response of our model to in¬ 
termittent driving allowed us to reproduce qualitatively 
an important aspect of the aging of contacts, namely 
the increase of the static friction with time of contact 
at rest. Overall, we believe our model reproduces many 
well-known features of real tribology, and gives a well de¬ 
fined model on which many assumptions and predictions 
of phenomenological theories (like rate-and-state equa¬ 
tions) can be investigated in depth. 
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Appendix A: Detailed Derivation of the Equations 
of Motion 

We first study the one-dimensional case, as presented 
in Fig. 9. 

The interface is decomposed in blocks of mass m, la¬ 
belled i and moving along horizontal rails hi. The ac¬ 
tion of the dashpot is to resist the change in (j)i — hi 
via viscous friction, with a resulting force on hi given 
by Tjudt{4>i — hi). The blocks move in a medium with 







FIG. 9. Sketch of the one-dimensional viscoelastic inter¬ 
face model. The interface itself (bold black line) consists in 
blocks set at discrete lattice sites i, f + 1,... along the x axis 
(in higher dimensions, x G R'*). The blocks can evolve along 
the 2 axis: they are identified by their locations hi, hi+i,... 
(empty squares) and bound together via a combination of 
springs (feijfe) and dashpots {rju). The viscoelastic interac¬ 
tion introduces an internal degree of freedom (f>i, represented 
by full squares (blue). Driving is performed via springs ko 
linked to the position w = Vt (thin purple lines). The dis¬ 
order force (red) for the site i derives from a disordered 
energy potential Ef'^{z) (grey). 


some effective viscosity rj and we study the overdamped 
regime, mdfhi r]dthi. As each block is described by 
two degrees of freedom hi and (fi, the time evolution is 
governed by two equations. We now provide a pedestrian 
derivation of the equations, for the sake of completeness. 
The first equation comes from the force balance on hf. 

V9thi =ft''{hi) -b ko{w - hf) + ki{h^+i - hi) 

-f- ki)hi—i hi) -\- hi) -t- k2(,4^i—i hi) 

(Al) 

The second equation is derived from the force balance on 

0 = k2{hi+i - (fi) + r]^dt{hi - cfi) (A2) 

where we assume that the internal degree of freedom (fi 
has no mass. Similarly the force balance on (fi-i yields: 

0 = k2{hi - (/)i_i) -b r]udt{hi-i - (ft-i). (A3) 

In order to let the Laplacian term fc 2 (^i-i-i ~ 2/ii -b /i^-i) 
appear, we introduce the variable 

Ui = (fi-hi +hi-i-(fi-i, (A4) 

which represents the elongation of the dashpot elements 


connected to site i. We inject (Eq. A2) into (Eq. Al) to 
get rid of the time derivatives, and we subtract (Eq. A3) 
from (Eq. A2) to obtain (Eq. A6): 

■qdthi =ft%hi) -b ko{w - hi) 

-b {ki -b fc 2 )(/ii+i — 2 hi -b hi-i) — k 2 Ui (As) 
hudtUi =k2{hi+i - 2hi -b - k2Ui. (A6) 

A more elegant notation using the Laplacian operator 
is: 

Tldthi = /^“(/li) + ko{w - hi) -b kiWjhi + k2{Vjhi - m) 

rjudtUi = k2iy'lhi - Ui). (A7) 

To generalize this to higher dimensions (on a square 
lattice), one simply has to connect each block hi to its 
neighbors via viscoelastic elements, using a single orien¬ 
tation per direction. The equations obtained are exactly 
(Eq. A7) if we reinterpret the label i as referring to d- 
dimensional space, the Laplacian as the d-dimensional 
one, and the Ui variable as: 

d 2d 

Ui = 'y ~ hj) A y ( {hji — 4>ji), (AS) 

j'=d+l 

where indices j denote the d first neighbors, connected 
via a dashpot followed by the spring k 2 (and ki in paral¬ 
lel) and indices j' denote the last d neighbors, connected 
via the spring ^2 followed by a dashpot (and fci in par¬ 
allel). 

We study the mean field limit via the fully connected 
approximation. In practice, each block position hi in¬ 
teracts with the positions of all other blocks via — 1 
springs of elastic constant ki/N [N being the number 
of blocks in the system) and via — 1 viscoelastic ele¬ 
ments (i.e. spring in series with a dashpot). As usual for 
fully connected systems, the final equation for the site i is 
obtained by replacing any occurrence of Ahi with h—hi. 


Appendix B: The Numerical Methods 

In previous works [8, 35] we showed how to trans¬ 
late the stochastic dynamics into a master equation (or 
Fokker-Planck equation), but only in the case in which 
E —^ 0. Here we generalize the procedure to any velocity 
V. Our Cython code is fully available [36]. 

It is useful to describe the dynamics using the local 
variable. Si: 

5i = l- ko{w - hi) - (ki -b k2)ih - hi) -b k2Ui, (Bl) 

which represents the amount of additional stress that a 
site can hold before becoming unstable: as long as Si > 
0, Vf, the block is stable and its position, hi, remains fixed 
when = 0 the block jumps to the next pinning well at a 
random distance z: hi ^ hi-\- z. To define the dynamics 
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in terms of 6 variables, we need to split it in a fast part 

, and a relaxed one 6 ^: 

S[ = 1 — ko{w — hi) — {ki + k 2 ){h — hi) 

= k 2 U„ (B2) 

such that 6 i = 6 f+ 6 ^. The blocks dynamics is controlled 
by three processes: 

(i) The avalanches. When a block is unstable {Si < 0), 

it moves to the next pinning wells: 5f >—>■ Sf + 
z{ko + ki), with z drawn from g{z). Each jump z 
is followed by a stress redistribution 5^ ^ ~ 

zki/N. These drops can trigger other instabilities. 
The characteristic duration of such an avalanche is 
To = r]o/ uiax{ko, ki,k 2 ). 

(ii) The driving. External driving increases over time: 
for instance w = Vt in the uniform driving proto¬ 
col. In terms of <5, driving means that over a time 
step dt, Sf I—>■ Sf — koVdt. This driving can happen 
until a new instability is triggered (step (i)). The 
characteristic time scale of driving is td = z/V. 

(iii) The relaxation. In absence of instabilities, the hfs 
are constant, and Eq. 2 reduces to: 

= k2(h-hi) -T {S^{to) - k2(h-hi)) e~^^^ 

(B3) 

where to is the time at which the last avalanche 
occurred. Note that hi does not evolve during re¬ 
laxation or driving, so that relaxation can happen 
until a new instability is triggered (step (i)). The 
characteristic time scale of relaxation is = gu/k 2 - 

In this paper we study the case tq ^ td and focus 

on the case of the mean field model where fluctuations 
vanish and the description of the system via a simple 
probability distribution becomes exact. Indeed, the sole 
distribution P{S) does not provide enough information to 
fully characterize the system and its evolution. We con¬ 
sider the joint probability density distribution P{S^, S^). 
The quantity P{ 6 ^,S^,t)dS^dS^ represents the proba¬ 
bility for a site drawn at random to have a particular 
set of values of 6 ^, 6 ^ and can be computed numerically 
starting from the dynamical rules that apply to the (5’s. 

To be concrete, we discretize P{ 6 ^, 6 ^) with a bin e. 
The distribution probability is then a matrix Pij where 
P{S^ = £i,S^ = ej)dS^d 6 ^ = Pij. We use a time step 
dt = e/koV and define the constant k = k^ + ki + k 2 . 

The finite velocity is expressed through the fact that 
every time there is some driving of ic by a quantity 
dw = Vdt, there is also relaxation during a time dt. 
In particular, we define the relaxation factor i?(dt) = 
I — For the avalanches two crucial quanti¬ 

ties should be defined: (I) the fraction of unstable sites 
-funst = ^^i,j\i+j<0 ~ d (^) 

stress redistribution P^°l^z{ki + ^ 2 ) that follows the sta¬ 
bilization of the unstable sites. When P*°gjz(/ci-|-/c 2 ) > 1, 


the avalanche increases geometrically over the time steps, 
this is why it is practical to define a “critical” value 

pc 1 

0 z(ki-\-k 2 ) ' 

The sketch of the algorithm is the following. 

• Relaxation process: 

— Compute joo(*), the bin associated to the fully 
relaxed state, = ^(S^ — S^). 


joo(*) = Int k2 




— Relaxation corresponds to shift [? 
] Pij to Pi,j+shift (where Shift = 
Int [(joo(*) — j)^(dt)]), set r = 1 and 

perform the Avalanche process. 

• Avalanche process: consists of driving and jumps. 

— Driving: 

Pi,3 ^ Pi,3 + {Pi+i ,3 ~ Pi,3) ^ (B^I) 

— Compute T’unst(j) = '^i'\i'+j<:o Pi',j 
— Jumps: 

Pi,3 ^ Pi,3 + -9 ( Punst(j) (B5) 

K \ K ) 


Pij i — 0 if i H- j <C 0 

- Compute Po = Y.i=-j Pi ,3 

~ If To > Pq, set r = I and perform the 
Avalanche process again. 

— Else 

* Compute = Y.j Punstij) 

*If > PS/m set r = 

min(l, ) and perform the 

Avalanche process again. 

* Else, perform the Relaxation process. 

Here we used the exponential distribution with z = 0.2 
and an upper-length cutoff g{z) = l/ze“^/^0(z)0(lO — 
z), where 0 is the Heaviside function. We always used 

Vu = I- 

Above we consider that the stationary regime is 
reached when the relative variation between the last three 
values of a does not exceed 1%. Below we ask that 
the stress drops Act of three consecutive Global Shocks 
change less than 1%. Finally the critical value is found 
by decimation: starting with a very small Vmin = Ie“^ 
and very large Vmax = 1000, simulations are repeated be¬ 
tween the two boundaries until the relative difference be¬ 
tween those two boundaries is less than 1%. This means 
that a larger absolute tolerance is allowed for larger val¬ 
ues of EC 
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FIG. 10. (Color online) Convergence of the stress to its 
limit value <7{V), as the binning e decreases. We selected 
velocities V close to but larger than for the five values of 
ko studied (from top to bottom: decreasing ko and increasing 
V). Convergence is reached for the three largest values of ko. 


An issue limiting higher precision is the presence of 
numerical instabilities when approaching to the critical 
point. In Fig. B the dependence of tr on the binning e 
is shown. For large values of fcp (0.3, 0.1 and 0.03), it 
is clear that we converge towards a plateau <j{y) when 
the binning decreases. The small values of fco represent a 
challenge, because the precision required from the algo¬ 
rithm is ^ zfcp. However, comparing the expected values 


(dashed lines) and the behavior at larger fcg’s, it can be 
expected that smaller binnings would produce the ex¬ 
pected results. 


Appendix C: Analytical solution for the stationary 
regime 

Let us now solve the model of Fig. 2(b). We call /°(t), 
and /^(t) the forces coming from the branches with 
the fco) k 2 springs. The coordinate hi will jump 

io hi-\-~z each time 


+ + = (Cl) 

Between jumps, the different forces behave as 

f{t)=fa+kiVt (C2) 

f\t)=fl+koVt (C3) 

fit) = Vv + {fa-'^Vu)exp(^-^^ (C4) 


where the a sub-indexes on the right indicate values of 
the forces right after the jump. This expression hold 
until the next jump that occurs at time t = z/H. At this 
moment, the forces must satisfy Eq. (Cl), and we obtain 


fa + klZ 


(C5) 

/a + koZ 


(C6) 

Vil + ifl - Vr]u)exp 

k2Z \ 

. yvu) 

(C7) 


where the b sub-index stands for the values right before 
the jump. In addition, we have = —f^ since the 
average value of must vanish. Also, since the dash- 
pot is rigid at the jump, we get — fa = ^^ 2 - From 
all these equations all forces (/^°, /°, fa) can be 

calculated. In particular we are interested in the aver¬ 
age friction force cr which is given by cr = (/° -|- /°)/2. 
Through a straightforward elimination procedure we get 
the Eq.(4) given in the main text. 
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